if (! require(edgeR)) {
   source("https://bioconductor.org/biocLite.R")
   biocLite("edgeR")
   library(edgeR)
}

if (! require(DESeq2)) {
   source("https://bioconductor.org/biocLite.R")
   biocLite("DESeq2")
   library(DESeq2)
}

data = read.table("/home/cuilj_ypd/wd/22.03.18_meta_ypd/5.annotation/unigenens.count", header=T, row.names=1, com='')
col_ordering = c(7,8,9,1,2,3)
rnaseqMatrix = data[,col_ordering]
rnaseqMatrix = round(rnaseqMatrix)
rnaseqMatrix = rnaseqMatrix[rowSums(cpm(rnaseqMatrix) > 1) >= 2,]
conditions = data.frame(conditions=factor(c(rep("ck", 3), rep("maize", 3))))
rownames(conditions) = colnames(rnaseqMatrix)
ddsFullCountTable <- DESeqDataSetFromMatrix(
    countData = rnaseqMatrix,
    colData = conditions,
    design = ~ conditions)
dds = DESeq(ddsFullCountTable)
contrast=c("conditions","ck","maize")
res = results(dds, contrast)
baseMeanA <- rowMeans(counts(dds, normalized=TRUE)[,colData(dds)$conditions == "ck"])
baseMeanB <- rowMeans(counts(dds, normalized=TRUE)[,colData(dds)$conditions == "maize"])
res = cbind(baseMeanA, baseMeanB, as.data.frame(res))
res = cbind(sampleA="ck", sampleB="maize", as.data.frame(res))
res$padj[is.na(res$padj)]  <- 1
res = as.data.frame(res[order(res$pvalue),])
write.table(res, file='unigenens.count.ck_vs_maize.DESeq2.DE_results', sep='	', quote=FALSE)
write.table(rnaseqMatrix, file='unigenens.count.ck_vs_maize.DESeq2.count_matrix', sep='	', quote=FALSE)
source("/home/cuilj_ypd/wd/22.02.26_meta_test/P2.Taxonomic_profiling/script/R/rnaseq_plot_funcs.R")
pdf("unigenens.count.ck_vs_maize.DESeq2.DE_results.MA_n_Volcano.pdf")
plot_MA_and_Volcano(rownames(res), log2(res$baseMean+1), res$log2FoldChange, res$padj)
dev.off()
